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Abstract 

An intcrfacial approximation of the streamer stage in the evolution 
of sparks and hghtning can be formulated as a Laplacian growth model 
regularized by a 'kinetic undercooling' boundary condition. Using this 
model we study both the linearized and the full nonlinear evolution 
of small perturbations of a uniformly translating circle. Within the 
linear approximation analytical and numerical results show that per- 
turbations are advected to the back of the circle, where they decay. 
An initially analytic interface stays analytic for all finite times, but 
singularities from outside the physical region approach the interface 
for t —^ oo, which results in some anomalous relaxation at the back 
of the circle. For the nonlinear evolution numerical results indicate 
that the circle is the asymptotic attractor for small perturbations, but 
larger perturbations may lead to branching. We also present results 
for more general initial shapes, which demonstrate that regularization 
by kinetic undercooling cannot guarantee smooth interfaces globally in 
time. 

PACS: 47.54.-r 
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1 Introduction 

Propagating fronts in Laplacian growth occur naturally in quite a number of 
physical problems including viscous fingering [H El El HJ [5] , electro-chemical 
growth, dendritic crystal growth for small undercooling [6l[71[8], and void 
migration in a conductor [9l \T0\ lllj . More recently, it has been shown 
that this class of problems includes the 'streamer' stage of electric break- 
down [HmSlEKliKinKIZKISKIS], which win be described below. A central 
issue in these problems is the stability of curved fronts. In a limiting case, 
most of these models reduce to the classic Saffman- Taylor problem [T] , which 
is known to be ill-posed \20\ I21j . Numerical as well as formal asymptotic 
results m El El E2] suggest that one branch of steadily propagating finger or 
bubble solutions in a Hele-Shaw cell is stabilized by surface tension regular- 
ization, though only recently some mathematically rigorous results [231 
are available to justify nonlinear stability to small disturbances in the special 
case of a nearly circular bubble. Besides surface tension, other regulariza- 
tions [HI [ini ttH [25] have also been analyzed. In the present paper we study 
both the linear and the nonlinear initial value problem for one such regular- 
ization, in particular, the stability of a steadily propagating circular shape. 
This regularization is called kinetic undercooling in the crystal growth con- 
tex10, but has a different physical interpretation for streamers. 

During the streamer stage of electric breakdown the discharge paves 
its way through a nonconducting medium, leaving behind a weakly ionized 
conducting channel. The basic growth mechanism is impact ionization due 
to electrons strongly accelerated in the local electric field. In a sufficiently 
strong field, a thin space charge layer forms around the head of the streamer. 
This layer screens the field in the inner ionized region to a very low level, and 
the growth of the streamer is driven by the electrons moving and multiplying 
in the strong self-enhanced field ahead of the curved ionization front. 

For sufficiently strong external fields, the thickness i of the electron layer 
is small compared to the radius R of the streamer head. Therefore Lozansky 
and Firsov suggested (mainly in the Russian literature, but also in that 
this layer can be modeled as an interface separating the ionized from the 
non-ionized region. Probably, the idea is even older, since a similar concept 
was already proposed by Sammer in the German literature in 1933 |27| . 
However, a deeper study of the implications of this concept started only 

^In crystal growth, kinetic undercooling is likely to be the more important regulariza- 
tion compared to the Gibbs- Thompson effect for large undercooling where the Laplacian 
growth model becomes questionable. This happens when the time scale over which the 
interface evolves becomes comparable to the time scale of heat diffusion. 
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later [HI [HI [H Ell EQ] where the problem is placed in the context 
of other Laplacian moving boundary problems. The validity of the moving 
boundary approximation for negative streamers is discussed in [15] for simple 
gases like pure nitrogen or argon, and in [31] for air. The dimensional 
analysis and the proposed regularization mechanism of the moving boundary 
problem for negative streamer ionization fronts are discussed in detail in our 
previous papers [IHl dZl dSl [TO] . 

In dimensionless form, the model is defined as follows. The normal veloc- 
ity Vn of the interface is given by the drift velocity v of the electrons, which 
is proportional to the local electro-static field E = — In appropriate 
units it takes the dimensionless form: 

Vn = -n-B+, (1) 

where the super-script denotes the limiting value as the interface is ap- 
proached from the exterior (the non-ionized region) and n is the outward 
normal on the interface. Outside the streamer the electric potential 93 obeys 
the Laplace equation: 

Atp = 0. (2) 

An analytical and numerical analysis of the underlying physical model for- 
mulated in terms of partial differential equations for charge densities and 
field suggests the interfacial condition 

^+ = en.(V(/.)+, (3) 

where 

Far from the streamer, the electric field tends to a constaniH 

E = -Vv3^ - x + o(l/|x|), (5) 

where x is the unit vector in x-direction. Eqs. ([I|)-([5|) define our model. 

In two dimensions, a simple solution to the free boundary problem posed 
by this model takes the form of a uniformly translating circle. Our previous 
work in [n\ \TE[ [T9] and the present paper are primarily concerned with 
the linear and nonlinear stability of this solution to small perturbations. 
It is to be noted that the circular shape differs from an actual streamer 
shape. However, the front half of a circle roughly resembles the shape of 
the front part of a streamer. Since growth of disturbances is found to be 
most pronounced in this advancing part of the interface, we expect stability 



^ A correction of order 0(l/|x|) to the electric field E can occur only if the streamer 
carries a net electric charge. We here concentrate on the analysis of streamers that are 
globally electrically neutral. 
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features found here to be qualitatively relevant for an actual streamer and 
more generally for curved fronts. 

In the special case e = 1, the linearized evolution of small perturbations 
can be determined exactly in our model [171[18]. The case of general e > is 
treated in part I [19] of this series of papers and in the present manuscript. 
In [19], we discussed the spectrum of the linear operator which results from 
the linear stability analysis of the circular solution. Restricting ourselves to 
an appropriate space of analytic perturbations we found a pure point spec- 
trum. Asymptotically in time, except for the trivial translation mode, all 
eigenmodes were found to decay exponentially in time. These eigenmodes 
are singular at the back of the bubble; nonetheless, as evidenced in the 
present paper, this singularity is not reflected by the actual linear evolution 
near the bubble back. The usual asymptotic form of the solution for large 
time: XIa^'^*/^^) where P\ is the eigenfunction corresponding to the eigen- 
value A, fails in a neighborhood of the rear of the bubble, though it holds 
elsewhere. 

In the present paper, we consider the initial value problem. For the 
linearized evolution, analytical results are obtained in the limit of strongly 
localized disturbances of the circle. Also the large time behavior of general 
perturbations can be studied analytically. Numerical calculations confirm 
these results. Together with the eigenvalue analysis of the first paper |19j . 
clear evidence of linear stability is presented. The full nonlinear evolution 
of a perturbed interface is calculated numerically. Our results suggest that, 
similar to linear evolution, small enough perturbations of a circular bubble 
grow in the front part of the bubble, but eventually decay as interfacial dis- 
tortions advect to the bubble rear. Nonetheless, when e is small but nonzero, 
the large transients in the linear regime make nonlinearity important even 
when the initial perturbation is exponentially small in e. Furthermore, when 
the perturbations are larger, the circle is no longer an attractor of the dy- 
namics and the propagating structure branches. For general initial shapes, 
we give some numerical evidence that the undercooling regularization condi- 
tion can not guarantee a smooth interface globally in time. For some initial 
conditions, the interface tends to develop a sharp corner in the back. Other 
initial conditions lead to the separation of the moving body into two parts. 

This paper is organized as follows. In section [2] we present equations 
derived earlier in a conformal map setting. Section [3] is devoted to the 
linear evolution of perturbations of the circle. Subsection I3AI recalls previ- 
ous results, and in subsection I3BI we present rigorous results on the growth 
of a strongly localized perturbation. We continue the discussion of local- 
ized perturbations in subsection I3CI and explain at an intuitive level how 
strongly localized perturbations are generically advected to the rear of the 
circle, increasing in amplitude in the front-half before decreasing in the back 
half. Mathematically, the advection is described by a one-parameter family 
of conformal maps which is a subgroup of the automorphisms of the unit 
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disk. The important role of this subgroup has been previously established 
for the exactly solvable case e = 1 [n\ [T8] . In subsection I3DI we discuss the 
anomalous behavior found at the back of the circle in the large time limit. 
In subsection I3EI we give arguments indicating that an initially analytic in- 
terface stays analytic for all finite times, but singularities initially outside 
the physical region of interest approach the back of the circle for t ^ oo. 
Provided the perturbation for t — > cxd stays analytic in the closed unit disk, 
except for the point —1, we in subsection I3FI prove that it asymptotically 
reduces to a constant. This implies that the perturbation just leads to a 
shift in space with respect to the unperturbed propagating circle. In sub- 
section [SHl we present numerical solutions of the linear evolution equations. 
These calculations support the asymptotic results derived in the previous 
subsections. For disturbances, not necessarily localized, we present evidence 
that on any part of the interface not containing a neighborhood of the bub- 
ble rear, the decay rate of the disturbance matches what is expected from 
the prior spectral analysis [19]. 

Section H] presents a numerical study of the nonlinear evolution for dif- 
ferent perturbations. We first consider perturbations of a circular bubble. 
It is shown that the circular bubble can be nonlinearly stable if the pertur- 
bation is small. However, when the perturbation is large enough, the front 
may start to branch. Furthermore, we study the nonlinear evolution for 
more general initial configurations. It is shown that the formation of a cusp 
precisely on the back side of the moving body can not be excluded. We also 
observe that the body might split into two parts. 

2 Equations resulting from conformal mapping 

As already explained repeatedly [El HTJ [THl [19] , we assume the streamer 
to be a simply connected compact domain V in the {x, y)-plane. The area 
of T> is conserved under the dynamics and equals vr in dimensionless units. 
Identifying the {x, y)-plane with the closed complex plane z = x + iy, we 
introduce a conformal map f{uj,t) that maps the unit disk U^j in the w-plane 
to the complement of T> in the z-plane 



The Laplace equation and the boundary condition (|5]) are incorporated 
in the definition of a complex potential ^{lo, t). 



z = f{uJ,t) 



a-i{t) 



a_i(t) > 0. 



(6) 




a-i(t) 



(7) 
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The remaining boundary conditions ([T]), ([3]) take the form 



Re 



dtf 



|9c./|Re $ 



Re 



UJ 



.\d^f\\ 
-eRe[ujd^<^], uj £ dU^ 



(9) 



(10) 



The problem reduces to solving these two equations, respecting the analyt- 
icity properties of / and <I>. 

A simple solution corresponding to a steadily translating circle is given 

by 



1 2t 

u 1 + e 



1 

UJ 



1 



(11) 



1 + e 



UJ . 



In physical space it describes a unit circle moving with constant velocity 
2/(1 + e) in x-direction. For small and smooth distortions of this circle, it 
is appropriate to look for solutions of the form 



/(w,t) = f(^){uj,t) + ri (3{uj,t) 



(12) 



where /3(u;, t) and x(^5 1) ^-re analytic in li^^ and r/ is a small parameter. Since 
the area is conserved, it can be shown that the residue 1 of the pole in (jlip 
remains unchanged to first order in rj. Substituting (fT2|) into equations 
(jlOP we in first order in rj find a system of two partial differential equations, 
from which x can be eliminated. The final equation for (3 takes the form 



-Ce = - <9<^ (w^ 



1) ujd^ + ed^^ujdr + dr - , 



where we introduced the rescaled time variable 

2 



1 + e 



t. 



(13) 
(14) 

(15) 



Eqs. ()13|) . ()14p determine the linearized evolution that will be discussed in 
Section [3l We will assume that the initial interface is analytic, i.e., that all 
singularities of /3(a;,0) are outside the closed unit disk U^j-, though much of 
the analysis is valid for a sufficiently smooth interface as weli. 



^Analyticity is not crucial, except in 
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3 Analysis of infinitesimal perturbations 



3A Summary of previous results 



In part I [T9j we have analyzed the eigenvalue problem, resulting from 
Eqs. p^ . via the ansatz /3(a;,r) = e^'^Px{io). We have shown that 
the spectrum is purely discrete and that the real part of all eigenvalues A„ is 
negative, except for the trivial value Aq = 0, which corresponds to a simple 
shift of the circle. An infinite set of real negative eigenvalues was found. All 
eigenf unctions, except for I3\^{uj) = const, are singular at w = — 1 at the 
back of the circle. Thus the expansion of a regular initial condition in terms 
of eigenfunctions has to break down in the neighborhood of w = — 1, which 
indicates that in that neighborhood some anomalous relaxation shows up. 
Furthermore, we found that as e j 0, any eigenvalue tends to zero and 
the corresponding eigenvector (3\{oo) tends to a constant. A similar behav- 
ior of the spectrum was found for a steadily moving circle in a Hele-Shaw 
cell with surface tension regularization [22] and this degeneracy is not unex- 
pected since the unregularized problem (e = 0) is mathematically ill-posed 
[201 E]. 

Here, we consider the initial value problem defined by Eqs. ([TH]) . 
Our analysis is guided by previous results [13 dH] on the special case e = 1 
where the general time dependent solution is known analytically; it is 



where the function is given by the initial condition, 

G(a;) = (2 + L^a^)/3(c^,0), (17) 

and r(r) is defined as 

r(r)=tanh^. (18) 

The properties of these solutions are discussed and visualized in detail in [TTJ 
[TH] . Here we in particular note that the essential time dependence of /3(a;, r) 
is contained in the transformation 

1 + ujT{t) 

(^{uj,T), < T < 1, defines a one-parameter family of automorphisms of the 
unit disk, with fixed points uj = ±1. The point w = 1 is stable, whereas 
cj = — 1 is unstable in the following sense: as r — > oo, i.e. T ^ 1, all 
the complex w-plane, except for a; = —1, is mapped into a neighborhood 
of C = +1- This results in an advective dynamics. Any perturbation not 
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centered precisely at a; = 1 is advected towards to = —1, where it vanishes 
asymptoticaUy. As r ^ oo, only a shift of the circle is left: 

hm /?(u;,t) = ^. (20) 

T— >00 Z 

However, it is to be noted that the limit is not uniform, and no matter how 
large r is, there is a neighborhood of w = —1, where P{u>,t) may change 
dramatically. We note that advection of distortions from the front to the 
sides has been observed in viscous fingering and crystal growth models with 
surface tension and has been derived from somewhat heuristically simplified 
models [H [32] . We further note that in the limit e ^ oo a purely advective 
dynamics results |18j : 

/3(w,r) =/3(C(u;,T(r)), e = oo . 

Expecting the automorphism (^{uj,T) and the resulting advective dynamics 
to play an important role also for e 7^ 1 we transform the PDE (|13p . ()14p 
from variables (w,t) to variables {(^,T), introducing the notation 

l3{L0,T)=mu;,T{T)),TiT)). (21) 

This results in the normal form of a hyperbolic PDE: 

ehiC, T)dT + 9c + (1 + e)dT } PiC, T) = , (22) 

where 

KC, T) = ~ = a;(C, T) [d^ a;(C, T)]-^ . (23) 



3B Localized perturbations; rigorous results 

Consider for general e > an initial perturbation that is centered at C = 
Q = e*'/'= for 7^ and has 'width' 7 in the sense that /3(e*'^'^+*'^^, 0) 
decays rapidly with |x| when |x| ^ 1. The decay rate will be specified more 
precisely below Eq. ([30]) . To study this problem, we first write (|22]) as an 
integral equation: 

^c(C,r) = c^A/.-iA(c,r)/3c(C,o) - £ Ps{C,s)h-'+'/^iCs)ds, 

(24) 

where /J^ and (3s denote derivatives of (3. Then integration by parts in s 
replaces Ps by (3 which is written as 

PiC,T)= (3^{c',T)dC' + fe,r). 

"'Co 
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Here Co = is a reference point in the tail of the perturbation chosen 
such that < ipo < ipc 1^ T^- We assume that ^""^^ is so large that (3{(^q,0) 
is negligible. 

Then, after some algebraic manipulation, we are able to rewrite (|24p as 
the following equation for 

G(x,r)^;3c(e''^^+'^^r). (25) 



Gix,T) = G'W(x,T)+ r r , ICi{x,x',T,s)Gix',s)dx'ds 

Jo J-t£^ 



+ 



V'e-i/'O 



IC2ix,x',T)Gix',T)dx' ^ G^'\x,T)+c\g] (x,T), (26) 



where, with the understanding that C = 6^'^'=+*'^^, C' = 6*'^'=+*^'^', 

nC'ii - e2)(i - r2)iA /(c - s)(i - sC)\-'+'/' 



K^iix, x',T,s) 



{(^_s)(l_^,)}iA V (1-^2) 



2s 



S-C SC-I S2 - 1 



/C2(X,X',T) = -ijC 



, (l + e)(l-T2) 

e(C-r)(i-rc) ' 



(27) 
(28) 



and 



G(°nx,T) = /.-V-(C,T)cV-/3c(C,0)- + (l + e)/3(C,0)C^/- 



e/i(C,r) 



eC/iV^(C,r) 



+ 



(1 - e^) 



/i"2+iA(C,5) /it(C,5) /3(Co,5)d^ 



(29) 



With G^^\x^T) considered knowr^ we determine the solution G{x,T) to 



the integral equation (|26[) for x G 
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, r G [0, To], where Xi? and 



To < 1 are some suitably chosen positive values independent of 7. Now it 
is clear from the expression for /Ci and /C2 that they are uniformly small in 
the 11. II 00 norm when is sufhciently small. We now choose the norm 



|G||^ sup sup l^(x)|G(x,T)| , 



(30) 



^Since f3{(o,T) cannot be determined without considering the full non-local problem 
on |i^| = 1, part of the expression (|29|l for G*-"' is not known. Nonetheless, if a disturbance 
is localized, the contribution to G''^' from l3{C,o,T) will be relatively small. In any case, 
in order to study the evolution in the x-scale, we are not prevented from considering G'°' 
as known. 
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where the positive weight function Vl^(x) obeys 

/Xr 
W-\x')dx' < C < oo x<XR- 
-oo 

For example, = for x ^ ^ 1 for x > would suffice for our 

analysis. We define G to be localized if HGH is finite, and Co) can be 
chosen such that P{Co,T) is negligibly small for T € [0, Tq]. Now it is clear 
from (j26p that the linear operator C has the contractive property 

||/:[Gi-G2]|| <C^||Gi-G2||. (31) 

It follows that there exists a unique solution to the integral equation (|26|) if 
is small enough and that for <C 1 

G(x,T)~GW(x,r), 

provided x and T are in the above specified range. For a perturbation 
localized in the sense given above, our result reduces to 

Gix,T) ~ h-^/' (e^'^o,r) (^e^^'^o+^^^o) . (32) 

We note that £.[G]{x,T) in general will not vanish for x ~^ co. This is the 
reason for restricting x to the interval given above and indicates that for 
T > the localized perturbation will sit on top of a dynamically generated 
delocalized background of amplitude ~ 

A detailed discussion of the result (j32|) will be presented in the next 
subsection. 



3C Localized perturbations; formal intuitive arguments 

It is useful to obtain the result (I32p through a more formal, yet intuitive, 
reasoning. This will also be helpful in our subsequent treatment of the long- 
time asymptotics in the anomalous region near the back of the bubble. We 
again restrict the analysis to the unit circle to = e*",a € M, or correspond- 
ingly to C = e^^jip ^ According to Eq. ()19p . the two angular coordinates 
a and tp are related through 

(l-r2)sinV' , , 

a = arctan — ;r- ■ — . (33) 

(1 -hT2)cosV - 2r ^ ^ 

Initially, (at T = 0), a and V obviously are identical. In terms of ^, the 
PDF (ESI) takes the form 



e hii;, T)dT + ^^f^^^ + i(l + e) St i /3 [e'^ ,t) = , (34) 
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where 

h{i;,T) = [d^ a) = ' ^ ' 

We now search for a solution that during its evolution stays localized near 
a fixed angle ■i/'c) with an angular width 7 ^ vr. We use the ansatz 

~P {e'^ = kocix.T) , (36) 

where again 

X = , 37 

7 

and (3ioc{XiT) is assumed to vanish rapidly for \x\ > 1- With this ansatz, 
Eq. (j34p takes the form 

e h{^, + 7x, T) 5^ + (dr H^c + IX. T)) + ^7(1 + e) dr 

■ Piocix,T) = . (38) 

For 7 <C vr we neglect the term 27(1 + e)c?T and the ^-dependence in the 
argument of h to find an approximate solution of the form 

PiocU, T) = /i-^/^(Vc, T) hocix, 0), (39) 



which is the same as (I32p . 

Before we evaluate this result we briefly discuss its limitations, as re- 
sulting from the present derivation. In view of the assumptions 7 <C vr, and 
Ixl < 1, the use of the zero order result h{tp) ^ h{ijjc) is justified provided 

(1 - T)^ + 4T sin^ y > 2T7X sin V'c + T{-yxf cos V'c • (40) 

This is valid for all times provided \ipc\ ^ 7, i.e., for initial conditions 
Plocix. 0) which essentially vanish in the forward direction ■0 = 0. For ■0c ~ 
the condition ()40p is violated if (1 — T) becomes of the order 7, and therefore 
the approximation becomes invalid in the large-time limit T(r) 1. This 
special role of perturbations in the forward direction is not unexpected since 
for such perturbations advection is ineffective. 

Neglecting the term ~ i7(l -|- e)^^ has more serious consequences. Sub- 
stituting into Eq. (j38p an ansatz of the form 



PUx, T) = (x, T) + 7 (x, T) + 0(7^) 

one finds that the result for [3^^^ violates the condition [3^^^ (x, T) ~ for 
Ixl ^ 1- A localized initial condition dynamically generates a delocalized 
contribution, with an amplitude proportional to 7/e^, in full accord with 
the rigorous discussion of the previous subsection. Again this result is not 
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unexpected since the eigenfunctions of the operator Eq. (fill) , are delo- 
cahzed. Assuming that we can expand an initially localized perturbation in 
terms of eigenfunctions we must expect that the balance of the expansion 
coefficients an e^" which for r = leads to localization, is destroyed by the 
time evolution. With these limitations in mind, we now discuss the result 



According to Eq. (fM|) . if expressed in the variable ( = exp {i{ipc + 7X)) 
the evolution of the perturbation is most simple. Neither its position ipc nor 
its shape PiociX}^) change. Only the overall amplitude h~^/^ varies with 
time. For < |^cj < 7r/2, i.e., if ipc is at the front half of the circle, h~^^'' 
increases up to a time Tm given by 

{T\Tm) + 1) COS Vc - 2T{Tm) = , (41) 

and then decreases again. For l-i/^d > vr/2,ft,^^/'^ decreases monotonically. 
For any -i/jc 7^ 0, we find the asymptotic behavior 

h-'/^^l^,,T{T)) ^ for r^oo. (42) 

sm ipcl 2 

For a perturbation centered precisely at the back of the circle (^/^c = tt), 
exponential relaxation 

h-^l^ {TT, r(r)) = (^i^) = e-^A (43) 

holds for all r. We recall that the localized approximation must break 
down it h~^^'' becomes of the order of the amplitude of the delocalized 
background. Nevertheless we will argue in subsection l3Dl that a contribution 
with asymptotic time behavior e~'^^^ generally shows up. 

Using Eq. ([33]) to transform back to u; = we see that the center 
adT^r)) is convected along the circle, reaching ztyr for r 00. A little 
calculation yields the velocity of this advection 

|:a,(r(r)) = sin ae(T(r)) . (44) 

This result has a simple interpretation. Recalling that we are working in a 
frame moving with the velocity v = x of the unperturbed circle, we identify 
the velocity (|44|) as the projection of v onto the tangent to the circle at the 
instantaneous location of the perturbation. 

In terms of Oc the overall amplitude of the perturbation takes the simple 
form 

^V.= f """-(^^;^^ V^\ (45) 
V sinac(O) J ^ ' 
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Figure 1: Evolution of a strongly localized perturbation for e = 1. Curves 
1, 2, 3 correspond to times r = 0, 1.84 and 4.59, respectively, a) Re (3 as 
function of the comoving angular coordinate X = (V' ~ V'c)/7- Broken lines: 
local approximation. Full lines: exact result with the shift of the circle 
subtracted, b) Perturbed interface in the physical plane in the system of 
local tangential and normal coordinates as explained in the text. 



It increases as long as the perturbation is on the front half of the circle and 
decreases on the backside. The maximum, reached for ac{T{T)) = ±7r/2, 
strongly depends on the initial position ac(0) = ipc- 

Defining the scale factor of the width of the perturbation as 



r = — 

dtp 



(46) 



we find 

r = &-'(fer(r)) = ^'°"''^';» . (47) 

Thus the width behaves similarly to the amplitude, except that for e ^ 1 it 
varies much less. For r — > oo it vanishes like e""^. 

So far we considered perturbations localized away from the tip V'c = = 
ac(0) of the circle. For ipc = 0, Eq. (I39p still holds for times such that 

1 - r(r) » 7 , 

cf. Eq. (j40p . It describes the initial increase and broadening of the pertur- 
bation. Advection, of course, is absent. For 1 — T ~ 7 the width becomes 
of order 1 and the local approximation clearly becomes invalid. 

On the qualitative level these results are most similar to the exact results 
found for e = 1 [17^118] and resemble the dynamics of a localized perturbation 
found in the context of viscous fingering [32] . 

The quantitative performance of the local approximation is illustrated 
in Fig. [H where for e = 1 the exact evolution of a localized perturbation is 
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compared to our approximation. From the exact result (I16p the contribution 
G{T)/2 representing a simple shift of the circle, has been subtracted. The 
initial condition is chosen as 

^fe*'^ 0) = - 

(e#_(l+^) e#c)2 

with 

^ =-- = — 
'^'^ 10 ' ^ 200 ■ 

Fig. [1^ shows Re/3 as function of x = (V' ~ V'c)/7 for three different times. 
Curve 1 shows the initial condition, where by construction the exact form 
and the approximation coincide. Curve 2 shows the perturbation when it is 
largest, in w-space being located near u; = —i. Curve 3 is taken at some later 
time. Evidently in this example the local approximation, (broken lines), is 
quite accurate. Very similar results are found for Im /3, which therefore is 
not shown. Fig. [TJ) shows the effect of this perturbation in physical space. 
In evaluating z = f{u),t), Eq. p2|) . we choose the amplitude t] = 0.007e~*'^'=. 
To combine the three curves into one plot, we introduced a time-dependent 
rotation of the coordinate system such that y' or x' are measured along 
the normal or the tangent to the unperturbed circle at the center of the 
perturbation, (i.e., at angle —ac (T(r)), since the inversion contained in the 
conformal map induces a sign change of the angles). In this representation 
the exact solution and the approximation cannot be distinguished within 
the resolution of the plot. We note that in physical space the shape of the 
perturbation varies due to interference with the unperturbed circle. 

3D Asymptotic relaxation near oj = —1 

In discussing the asymptotic relaxation we prefer to rewrite (I22p in terms of 
r, usmg (1 - T'^)dT = 2dr. Inserting the exphcit form ([23]) of h{C,T) and 
multiplying by (1 - T^) 2/2 we find 

e(C - T)(l - TC)dr + (2TC - ^(1 + r2)(l + C2))ac 

+ {l + e){l-T^)dr]p(,C,T) = 0. (48) 

Here T stands for 

T = T{t) = tanhr/2 = 1 - le'^ + O (e'^^) , 



cf. Eq. ([18]) . Keeping only the leading r-dependence in the coefficients of 
the derivatives, we reduce Eq. (fl5|) to 



{edr + m-Q? 9c + 4(l + e) e-^d^ /3(C,T(r)) = 0. (49) 
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For T ^ 1 we neglect the term e "^dr f3 to find 

/3(C,T(r))~e-^/^/3o(C) + 7o, (50) 

where /3o(C) and 70 depend on the initial condition /^(C, 0) and of course 
cannot be fixed by this asymptotic argument. 

Since the derivative 9^ in Eq. (H9]) is multiplied by (1 — C)^, the neglect 
of the term involving e~'^dr can be justified only for 7^ +1. In terms of 

= ^ = -1 + 2 + O (e-2^) (51) 

this implies that we deal with a neighborhood of uj = —1 that is contracted 
to this point like e~'^ . This range of uj is complementary to the region 
where an expansion in terms of eigenfunctions can be expected to be valid 
asymptotically. 

In the result ()50p the ^-dependence is suppressed by a factor e~^/^, which 
for e ^ 1, r — > 00, vanishes much faster than e""^. Thus ^-dependent 
corrections of order e~'^ will dominate the asymptotic relaxation at the back 
of the circle. Noting the presence of e~'^ in the coefficients of the differential 
equation, it is natural to determine the structure of these terms with the 
ansatz 

^(C,r(r)) = e-^/^/3(C,r)+7(C,r), (52) 

where 

00 

7(C,r) = J]7fc(C)e-^-. (53) 

k=0 

From Eq. ([l8|) with f3{(^,T{T)) replaced by ^{C,t) we find 

1 + e 4ci \ 



7(C, r) = CO + ci e-" + [c2 - 3-3^ J e"^" (54) 

1 + e 8C2 , (l + e)2 16ci ^ ^.3, , 4, 



+ 1^^-13^ W + (l- 2e)(l- 3.) (wFJ ^ 

where the Ck are integration constants. Generally 7a;(C) is found to be a 
polynomial in (1 — C)~^ of degree k — 1. In this analysis we assumed e 7^ 1/n. 
For e = 1/n the ansatz (I53p has to be modified. In particular a term 
proportional to re"""^ has to be included. We note that the exact result for 
e = 1 shows such a contribution [TS]. 

To transform our result back to w-space we introduce 

Slo = (l+w)e^, (55) 
j{5lo,t) = 7(c(M,r) -7(-l,r). (56) 
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Eq. ([IS]) yields 



1 + TTT '^'^ 1 1 

-- + ^6io + 0(e-") . (57) 



1-C 2-6uje-^ 2 4 
Thus 7(511;, r) has an expansion of the form 

00 

fc=i 

where the Ck again depend on the initial condition. For e ^ 1 the terms 
of order k < 1/e — 1 dominate over the contribution e~'^^''(3. For e < 1/2 
we therefore in a region of size |1 + lj| = O (e""^) near uj = —1 expect 
to see a very smooth asymptotic relaxation of the interface, with only a 
few coefficients depending on the initial condition. In contrast, for e > 1/2 
the asymptotic relaxation is determined by the term e~'^/'^/3o(C)! which will 
depend on the initial condition in a complicated way. For e = 1 this is 
illustrated in Fig. 5.2 of Ref. [18]. In the next subsection we will argue 
that the function /3o(C) picks up contributions due to singularities of the 
initial condition, which for r — > 00 are driven towards uj = —1. We finally 
note that the results discussed here resemble the behavior of the low order 
eigenfunctions (3\{uj). As shown in part I [19] of this series, these functions 
near u! = —1 develop a singularity of the form (1 + uj)^^^~^^, implying that 
the derivatives at u> = —1 exist for all orders k < 1/e + A. 

To illustrate our results we consider a perturbation centered at w = — 1. 
As initial condition we choose 

f3{C,0) = ^- , 7 = 1.05, (59) 
7 - C 

and we calculate the function 

B^,,),/^(-«-*.nr))-g(O.T(r))^ (60) 
f)(-l,r(T))-/3(0,T(T)) 

We expect to find the limiting behavior 

1 6*^^/2 1 

B{i;,T)r^- ^-1 = ^, e<-, (61) 

or 

I . ^ /3o(-e-^)-/3o(0) ^ 1 

B{'ip,T) T ^ o6 — , e>-, (62) 

/3o(-l)-/3o(0) 2 

respectively. Whereas /3o(C) depends on the initial condition, the limit (i6T]) 
is universal. The results shown in Figure [2] conform to these expectations. 
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a) b) c) 



Figure 2: Time evolution of the initial condition (j59p . a) \B{^lJ,T)\ as func- 
tion of ip/n for times r as given and for e = 0.1. b) Phase of B for the same 
values of r and e. c) Again \B(^l', r)| as function oiip/n, but now for e = 0.8. 



Fig. [2^ shows \B{iilj,T)\ for several values of r and for e = 0.1. It il- 
lustrates the approach to the limiting form l/cos(V'/2), which within the 
accuracy of the plot is in fact reached for r ~ 4. Fig. [2)3 shows the corre- 
sponding phase of B{iI),t). Here the approach to the limit is slower, but is 
definitely visible. Fig. [2t shows results for \B{ij),T)\^ e = 0.8. Here \B{iIj,t)\ 
seems to approach a limiting curve which clearly shows remainders of the ini- 
tial peak. (We should note that B{il), r) is symmetric: B{—iIj, r) = B*{tp, r), 
and that the peak at i/' = 0, of course, is rounded, which however is not vis- 
ible on the scale of the plot). We finally recall that the tp-iange shown here 
in terms of u; = corresponds to a small region near a = tt. Specifically 
for r = 4 it corresponds to vr < a < 1.08 tt. 

For the asymptotic relaxation our results predict 

[ 

/3(C,T(r))~<^ forr^oo. (63) 

This prediction is tested in FigureOby plotting results for ln[/3(— 1, T(r)) — 
/3(0,T(t)] as function of r for several values of e. The expected behavior is 
reasonably well observed. 

3E Analyticity of the interface 

If we assume the initial interface to be analytic, all singularities of the initial 
perturbation (3{lj,0) have to be outside the closed unit disk, We here 
argue that under the linearized dynamics the singularities stay outside U^j 
for all finite times r. For r — > oo they approach to = —1, and contribute to 
the anomalous e~'^^'' /9(C) behavior found in subsection I3DI 
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This argument is based on the recurrence relation for the coefficients 
6fc(r) in the Taylor expansion 



ijj 



(64) 



fc=0 



The evolution equation (fT3l) . ([H]) yields 
2 + e , 



2drbk 



1 + e ^' 

fc + 1 
1 + e + ek 



(65) 



(2 + e + eA;) bk+i - e{k - l)bk^i 



for A; > 1 . 



The singularities of (3{uj,t) are determined by the behavior of the ^^(t) in 
the limit k oo. For simplicity, we consider an initial condition (3{lo,0) 
whose singularity closest to |ti;| = 1 is a branch point at loq with behavior 
P{uj, 0) ~ const(tL'o — w)""^ for nonintegral a or a pole with a a non-positive 
integer. Then bk{0) for k ^ 1 behaves as 

bk{0) = const uj^'^k-'' (^1 + O Q 

We therefore make the ansatz 



bkir) = (-1)'= e 



-kf{r) 



9{r)k~ 



1 



1 



(66) 



where the factor (—1)^ is introduced since we expect the point u; = — 1 to 
play a special role. We will find that this ansatz is internally consistent 
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provided /c » fco, where ko increases with r. We shall conclude that for any 
r > 0, 

I e--^(^) I < 1, (67) 

if this condition is satisfied initially. This implies that the singularity re- 
mains outside the unit disk for all times and suggests that at least for initial 
conditions with a branch point, the interface will remain analytic. 

Substituting the ansatz (fBH|) into the recurrence relation we find 



-2drf + e 



-f _ J 



A 



1 



2drlng + {- + 1 



1 



a 



e-^ + e ^ 



(68) 



The leading order yields 



with the solution 



2drf 



-f _ J 



l-Ce 



where C is some integration constant. g(r) is determined by the next order: 



g{T) = ff(0)e(°-i-^)^fl-C2 



-2t 



(70) 



Checking higher orders in an expansion in powers of 1/k, one finds that 
neglecting such terms assumes k ^ ^ = ko. Combining our results we find 
the asymptotic behavior 



bkir) ~ (-1)* 



9(0) e 



a-l- 



l-Ce- 

Regularity of the initial condition enforces 



r 



k- 



(71) 



e-/(0) 




l + C 






1-C 



< 1, 



equivalent to Re C < 0. With the form (fB^ of /(r) this guarantees that 
condition ()67p . je"-^^^^! < 1, is fulfilled for all finite r. Thus for r < oo 
the singularities of (3{u),t) stay at some finite distance from the unit disk 
and the interface stays smooth, /(r) vanishes for r — > oo, indicating that a 
singularity reaches w = — 1. 

In the above ansatz (1660 . we assumed a particular type of branch point 
or a pole for /3(w, 0) as the nearest singularity. Multiple singularities of this 
type can be accommodated in this linear analysis using the superposition 
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principle. Other singularities can be accommodated as well by replacing 
by a more general k dependence. 
We now consider the limiting behavior of bk (t) for r — > oo more closely. 
Eq. dni) yields 

bkir) ~ (-l)^exp [2Cke-^] g{0) (A; 6"^)^+^"° fc^"^ . (72) 

This result, however, for r ^ oo is only valid for 

= A: e"^ > 1 , (73) 

i.e., for extremely large k. To extend the analysis to values rj = k = 0{1) 
we make the ansatz 

6fc(T)~(-l)H-7-i5(^,r), (74) 
which is motivated by Eq. (j72p . The recurrence relation (j65|) takes the form 

{2dr-2rid^) g{v,T) = (^k + O (^^^^ ~g {i] - r) 

k 



-2k e--d^~g{v,T) + o(j 



k 

= -2r?9^g(r?,r) + oQ^, 

or 

2d^~g{il,T) = o{^, (75) 

equivalently. Thus to leading order in l/Zc, g{'r],T) is independent of r and 
Eq. ([HT]) reduces to 

6fc(T)~(-l)^-fc-^i^o(fce--) . (76) 

Inspecting the terms of order 1/k one finds that this result asymptotically 
should be valid for e~'^ <C 1 and k > e'^/e^. 

The bkir), Eq. (f76l) . can be interpreted as coefficients of a Taylor expan- 
sion with respect to uj of the function /3o(C('-^) ^)) introduced in the previous 
subsection, cf. Eq. (I50p . To show this we again introduce 5uj 

to = —1 + 6uj e~'^ 

as defined in Eq. (I55p . and we approximately resum the Taylor expansion 
from k = r]Qe^ to infinity, using the result (j76p . 

oo , 

E hk{T){-l + 5u: e~y 
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oo 

~ J dk k~~''^gQ(k e~^) exp ^—k 5uj^ 

oo 

no 

This clearly is of the same form as the anomalous contribution in our pre- 
vious result (|5Up . The (unknown) function /3o(C) is given by the integral 
involving the (unknown) function goiil)- By construction the result (I76p is 
valid for large k and large r and therefore picks up the structure of the singu- 
larities for r ^ 1. We conclude that the anomalous contribution e~'^/'^/3o(C) 
is due to the singularities which approach w = — 1, as claimed above. 

We finally note that the leading singularity ~ (1 — uj)^^^'^^ of the eigen- 
function (3x{co) implies that the Taylor coefficients of e^'^f3x{Lo) for large k 
behave as 

M.) (I) V.(i + o(i)). 

where Ci is some constant. We thus recover the form (|76|) with go{r]) = r]~^. 
3F Rigorous analysis of the limit r — oo 

In the previous subsection we have argued that /3(u;,r) for r — > oo tends to 
a function f3oa{u)) that is analytic in any compact subset /C of lA^ \ {— 1}. 
Furthermore, the eigenvalue analysis [19] as well as the results of subsections 
I3BII3DI suggest that within the linearized theory a perturbation for r — > oo 
only leads to a constant shift of the circle. Assuming the existence of /3oo(w), 
this can be proven rigorously. 

We start from Eq. (fT3l) : C^fi = 0, rewritten as 

[(l-r2)5T- (l-w')9^] {l + e + eojd^)P{uj,T{T)) 

= (l-e) (l+a;2)5^/?(^,T(r)) , (77) 

where T = tanhT/2, (Eq. [TH]) . and we introduce the function 

G{uj,T) = {l + e + eujd^)(3{uj,T{T)) . (78) 

In terms of G{uj,T) the solution (3{lo,t) regular at a; = is given by 

Pico, t) = ^uj -1/^-1 J u '^/' G{lo ', T{T))duj ' (79) 



which generalizes Eq. (jl6p to e 7^ 1. We now write Eq. ()77p as 

[(1 _ t2) 5t - (1 - 0.2) d^] G(w, T) = H {uj, T) , (80) 
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where 



H{u;,T) 



1-e 



UJ 



1 + e 



x^/" G{xuj,T)dx 



(81) 



Noting that T) = C = + T)/(l + ojT) solves Eq. ([80]) for F = it is 
easily found that (j80|) is equivalent to the integral equation 



G(a;,r) = G(0,C) 
We now define 
Eq. ([82]) yields 

A(^,r) = G(o,c)-G(o,r) 



1 - w 



/2 



H[oj\ d^'- (82) 



A{uj,T) = G{uj,T) -G{0,T). (83) 



1-e f 1 l+a;'2 
^ 1 -a;'2 







A 



1 



1 + e 



xuj', — ^ 1 dx 

1 — U! C 



1 - W'C. 

dw', (84) 



where we have written out H explicitly. In view of the results of subsect [3El 
we now assume that lim Gioj.T) exists for w € /C. We further note that 

T-»l 

for r ^ 1 and uo ^ -1, both C and (C -uj'Q tend to 1. Eq. ^ 

reduces to the homogenous integral equation 



A(i^,l) 



'"7^1^ A(.',l)-i±i/xV.A(x.M)<<. 



duj' . 



(85) 



It is easily checked that for all e > the only solution of ()85p analytic in a 
neighborhood of u; = is the trivial one: 



Afw, 1) = 0. 



(86) 



To see this, we assume that the Taylor expansion of A(cj, 1) starts with a 
lowest order term w'^, A; > 1, 0^/0. Eq. ([85]) yields = 0, contradicting 
our assumption. 

We thus have shown that provided G(u;, 1) exists and is analytic for 
id G }C, the only solution to our problem is 



implying 



which for e = 1 reduces to Eq. (j20j) . 



G{uj,l) =G(0,1) 
G(0,1) 



(87) 



/3oo(w) 



1 + e ' 
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3G Numerical illustration 



In this section, we show numerical results of the linear evolution. We ap- 
proximately solve the PDE (I14p by truncating the series expansion ([6l 



fc=0 

at k = N. The ODE system for 6fc(r) has been given in Eq. (I65p 
^drbo = bi , 



2drk 



'k 



1 + e 

k + 1 
1 + e + ek 



(2 + e + ek) bk+i - e{k - 1)6^. 



for k>l. 



With the 6fc(0) given by the initial condition, the 6fc(r) can be determined 
recursively by the Runge-Kutta time stepping method. We choose the cut- 
off iV = 2000 in the simulation. Adaptive time steps are chosen which ensure 
that the difference between 4-th order and 5-th order Runge Kutta methods 
is within 10~^^. In the sequel we present results for 

5P{uj,t) = /3(c^,r)-/3(0,r). (89) 

The subtraction eliminates the overall shift of the evolving body. 

We first present results typical for a delocalized initial condition, choos- 
ing 

l3{uj,Q) =00^ . (90) 

Fig. m shows the evolution of Re oj5(3{lo,t) , uo = e*", with e = 1/10 or 

e = 1/2, respectively. In physical space, Re[a;(5/?] is the component of the 
perturbation normal to the unperturbed but shifted circle at angle —a. 
Panels a and b show that the qualitative behavior is quite similar for both 
values of e shown. Panels c and d give a detailed view on the state for several 
time steps; here an extended range of a is shown, so that the behavior both 
at a = and |a| = vr is clearly seen. For small times the perturbations 
increase in the front half |a| < 7r/2 of the circle and decrease in the back half. 
The maximum at a = increases and broadens strongly, whereas the other 
perturbations are shifted towards a = ibvr. At later times the perturbations 
decrease at the front half, while at the back a transient increase is observed 
which is due to the advection of the dynamically generated large amplitude 
of the perturbation towards a = vr. The results for e = 1/10 or e = 1/2 
essentially differ only in two respects. First, for e = 1/10 the perturbation 
at intermediate times is amplified much more than for e = 1/2. Second, 
for T = 2 and e = 1/2, remainders of individual maxima that initially 
are located at a 7^ 0, still can be seen near a = itTr (this structure was 
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(a) . 



-1 o/ti 
= 1/10 




(c) 



-1 -0.5 0.5 1 

Times 0, 0.5, 2, 6 for e = 1/10 



-1 

(b) 6 = 1/2 




(d) Same time steps for e = 1/2 



Figure 4: Evolution of Yle[uj 5P{uj,t)] for u = e*" for the initial condition 
(|9U|) . a) Overview plot for e = 1/10 and times < r < 11. b) The same for 
e = 1/2 and times < r < 11. c) Detailed data for time steps 0, 0.5, 2, 6, 
11 for e = 1/10; the angle is normalized as a/vr; orange areas are overlap 
regions plotted to make the structure at the back visible, d) The same for 
e = 1/2. 



very pronounced for e = 1 as discussed in [IH|)) whereas for e = 1/10 this 
structure is completely damped out and yields a broad maximum. 

Fig.[5]shows \5P\ as a function of a and r. It illustrates how the maximum 
of the absolute value of the perturbation is advected towards a = vr, where 
it decays. For e = 1/10 the behavior of \6(5\ is quite smooth, whereas for 
e = 1/2 some small scale structure is observed near a = ibvr. 

Outside a neighborhood of a = vr we expect to see asymptotically ex- 
ponential relaxation: 6(3 ~ e'^^'^. From the results given in paper I [TH] in 
Fig. m we expect Ai ^ -0.546 for e = 1/10, and Ai « -0.905 for e = 1/2, 
respectively. These predictions are tested in Fig. EJ Since the maximum of 
\Sl3\ advects along the circle, we plot In |5/3max(T)| as a function of r, where 

|5/?max(r)|= sup |(5/3(e^°,r)|. 

ae[7r+(5i ,7r — (5i] 
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-1 a/ji 

(a) € = 1/10 



-1 o/ti 

(b) e = 1/2 



Figure 5: Evolution of \5[3{uj,t)\ for the initial condition (I90p with (a) e = 
1/10 and (b) e = 1/2. 



We choose (5i = vr/S. For smaller values of it needs larger values of r to 
reach the asymptotic behavior. We fit the curve for data at 10 < r < 11. 
As Fig. [6] illustrates, the expected asymptotic behavior is observed. 
We now consider a more localized initial condition: 



7 



+ 

7 — u; 7' 



r 



UJ 



-y 

2 ^ 

j=0 



7/ V7 



(91) 



With the choice 7 = l.le*'^/^'^ it shows two fairly sharp peaks centered 
symmetrically close to a = 0. Similarly to Fig.Hl Fig.[7]shows Re(u;(5/3(ci;, r)) 
for e = 1/10 and e = 1/2. For e = 1/10 we rescaled the amplitude in panel 
c by a time dependent factor oo(t) in order to show all curves in the same 
plot. Panels a and b show that the time dependent shift of the structure 
is essentially independent of e, implying that advection is determined by 
the automorphism = (^(o;, r(T)). Panels c and d illustrate that also the 
detailed structure at given time r is fairly independent of e, but for e = 1/10 
the amplitude at intermediate times is enhanced much more than for e = 1/2 
(cf. the rescaling factors ao(T) given in the figure caption). Fig.[8]shows |(^/3| 
as function of a and r, similarly to Fig. [5l Again the advection of the 
maximum towards a = vr, its increase as long as it is in the front half, and 
its final decay in the back half are clearly seen. 

In summary, all numerical results presented here and in previous sub- 
sections support our analysis. 



4 Nonlinear Evolution 

In Sect. [3] we discussed the solution of the linearized evolution equation. 
Here we seek to determine the effect of the nonlinearity. In Subsect. I4AI we 
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slope = -0.54 slope = -0.90 




(a) e = 1/10 (b) e = 1/2 



Figure 6: logsupQ,gj_^^^^ \6P{a,T)\, 6i = vr/S as a function of r for the 
data presented in Fig. O a) e = 1/10, the line has slope Ai(l/10) = —0.54. 
b) e = 1/2, the hue has slope Ai(l/2) = -0.90. 



consider small perturbations of the circle. Subsect. I4BI presents examples of 
the evolution of more general initial shapes. 

To calculate the nonlinear evolution in a large range of time is difficult. 
Using a Fourier representation of the interface it for large r needs wave 
numbers of order e"^ S> 1 to resolve the collapsing region near uj = —1. 
Modes of large wave number can also be expected to play an important role 
at the front part of the bubble. Approximating a small region near lo = +1 
as planar, we may invoke well known results \28\ [29l [33] on the instability of 
a planar interface: in linear approximation the amplitude of a Fourier mode 
of wave number k increases like e'^^'^^^, where 

s{k) ^ 



1 + ek 



Thus with the present regularization all Fourier modes are unstable, whereas 
with curvature regularization only a finite unstable band exists. The strong 
increase of a perturbation localized near uj = +1, as discussed in Subsect. I3C[ 
can be considered to result from this instability of modes k ^ 1/e. Nonethe- 
less, despite stringent demands on resolution and time steps, we believe that 
the results presented here exemplify the nonlinear effects. The numerical 
methods used to solve the nonlinear equations ^ and (jlOp are summarized 
in the Appendix. We use a Fourier representation with cutoff kmax = 
and we solve the resulting system of ordinary differential equations with a 
4th order Runge-Kutta method with time step At. (A^ and At are given 
in the figure captions.) The numerics abruptly breaks down at some time 
imax(e). The time range shown in the figures therefore depends both on e 
and on the initial condition. 
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(c) e = 1/10 (d) e = 1/2 



Figure 7: Evolution of Yie[uj 5(3{u),t)] for the initial condition (j9ip with 
7 = l.le~^^^^'^ for e = 1/10 and e = 1/2. In panel c), Re{uj6P) is scaled by 
factors aoir), where qo(0) = 1, ao(l) = 0.03, ao(2) = 0.007, qo(4) = 0.025, 
ao(6) = 0.05. 



4A Small perturbations of the circle 

We here consider perturbations ij j3{uj, 0) of the circle, with ^ 1. In order 
to compare with the linear evolution we define the nonlinear counterpart to 
5P{oj,t) = P{uj,t) -/3(0,r) as 

/(a;,t)-/(0,t) 
dfini{(^,T) = , (92) 

V 

where t = ^-^t, Eq. (fTSl) . and f{uj, t) is defined in Eq. ([6]). For — > 0, 6Pni 
reduces to 5f3. 

We first consider the delocalized initial condition (|90p : f{uj,0) = rjoj^. 
Even for very small rj it is not obvious a priori that the nonlinearity is 
unimportant. As recalled above, perturbations at the front may increase 
dramatically, and the collapsing region at the back, where an eigenmode 
expansion is bound to fail, also might be quite sensitive to nonlinear effects. 
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(a) e = 1/10 



(b) 6 = 1/2 



Figure 8: The same evolution as in Fig. [TJ plotted is now |(5/3(a,r)| as a 
function of q and r. Left column: e = 1/10, right column: e = 1/2. Note 
the large difference in the scale of \6P\, reflecting the strong e dependence 
of the amplitude. 



We therefore in Fig. [U] show (5/3„;(+l,r) and l,r) for e = ^ or ^ and 

several values of r/. It is seen that for very small values of rj the nonlinear 
theory essentially reproduces the results of the linear approximation. Devi- 
ations outside some initial time range become visible for i] > lO""^, (e = jq), 
or 7] > 10~^, (e = respectively, but even then the shape of the curves 
is similar to the linear approximation. This suggests that also in the for- 
ward and backward regions the nonlinearity for small perturbations does 
not qualitatively change the results of the linear approximation. 

The results shown in Fig. [10] support this conclusion. We here plot 
i?e(e*"5/3„;(e*°, r) as function of a/vr, for values of r where a deviation 
from the linear approximation is visible. We observe that the nonlinearity 
essentially influences the amplitude but not the shift of the perturbation. 
The overall structure is most similar to the linear approximation. 

Such results are also found for other delocalized perturbations of type 
r/o;". Also more localized perturbations behave similar. For the initial con- 
dition dHI]), (3{uj, 0) = ^75^ + ^7:^ this is illustrated in Fig. dH We again 
observe that the nonlinearity essentially influences only the amplitude of 
the perturbation, but leaves the qualitative structure almost unchanged. 
All these results suggest that the circle is the asymptotic attractor for weak 
perturbations. 

For larger initial perturbations it is unlikely that the circle is recovered 
asymptotically. Rather we may observe branching. This is illustrated in 
Fig. [12] with the initial condition P{uj,0) = —0.03 to^. Panel (a) shows 
snapshots of the interface in physical space z = x + iy = f{io,t), with 
e = as resulting from the nonlinear evolution. For comparison panel 
(b) shows the linearized evolution, and panel (c) shows the result of the 
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Figure 9: Initial condition: /3„;(u;,0) = w^, plotted are (5/?(l,r) and 
6P{ — 1,t) as a function of r for e = 1/10 and e = 1/2 for different rj. 
{N = 256, At = 0.001) 
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(a) e = 1/10 



(b) e = 1/2 




(c) e = 1/10, r? = 10"^ (d) e = 1/2, r; = 10^^ 



Figure 11: Initial condition: /5(tj, 0) = + \ ^ ; plotted are 

supQ,g[_^+5^^7r-(5i] ^1 = as a function of r for (a)e = 1/10 

and (b)e = 1/2 for different r/. Evolution of Re[u; 5/3(u;, r)] for = e*" for 
(c)e = 1/10 and (d)e = 1/2 at different r. (iV = 256, At = 0.001) 
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(a) e = 1/10, nonlinear 



(b) £ = 1/10, linear 



1.5 



-1.5 




-1.5 -1 -0.5 0.5 1 1.5 2 2.5 



(c) 6 = 0, nonlinear 

Figure 12: Initial condition: P{uj,0) = —0.03 u^. (a) and (b) show the 
nonhnear and the hnear evolution respectively for e = 1/10. (c) shows the 
nonhnear evolution for e = 0. {N = 512, At = 0.00025) 
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E = 0.2 




unregularized model e = 0. Snapshots are taken at times t = 0.05n, where 
n = 0, 1, . . . , 12 in panels (a) and (b), and n = 0, 1, . . . , 5 in panel (c). Clearly 
the cusps which in the unregularized model occur for t ~ 0.25 for ^ = jq 
are suppressed both according to the linear and the nonlinear evolution. 
A qualitative effect of the nonlinearity is observed for t > 0.1. Whereas 
the linear approximation develops shoulders connected by some flat part of 
the interface, the nonlinear evolution results in two branches separated by 
a valley. Since the bottom of the valley moves slower than the tips of the 
branches, the valley is likely to evolve into a deep fjord. 

To close this subsection we briefly consider the range of validity of the 
linear approximation. As is evident from Fig. [U for a given initial condition 
this range strongly depends on e. The results of Subsect. I3CI suggest that 
it might decrease exponentially: ij < const g-™^***/"^^ where the constants 
might depend on the initial condition. To test this hypothesis, we for the 
initial condition /3{lo, 0) = 7]uj^^ compared the linearized and the nonlinear 
evolution for values r] = 0.001 • 2"", n = 0, 1, . . . , 5 and e in the range 0.12 < 
e < 0.24. We specifically calculated the absolute value of the difference 
dwPni ~ dul^ in forward direction a; = 1. We choose the derivative since it 
prominently shows up in the nonlinear equations @, (|10|) . The results for 
e = 0.2 are shown in Fig. [TST a). We observe that after some initial rise 
this difference saturates at some ?/-dependent plateau, where the plateau 
value strongly increases with -q. Eventually it decreases again, in agreement 
with the expectation that for the small perturbations rj, the circle is the 
asymptotic attractor. Interpolating among the plateau values we now for 
each e determined a value r/*(e) where the plateau value equals 0.02 at r = 4. 
Fig. 113( b) shows Inr/* as function of 1/e. As expected, it shows an essentially 
linear decrease. This supports the hypothesis that the range of validity of 
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the linear approximation, and presumably also the basin of attraction of the 
circle, decrease exponentially with increasing 1/e. 




Figure 14: /(w, 0) = - 2u; + 0.5u;2 (N = 512, At = 0.0005). The solutions 
are shown at t = 0.2n. 



4B Examples of the evolution of general shapes 

For general initial conditions the time evolution may lead to a breakdown 
of the model by two different mechanisms. First, a global breakdown oc- 
curs at time tc where the mapping looses the property of being one-to-one: 
fie^^^jtc) = f{e^°'^,tc), a\ ^ 02- Clearly for t > tc the model becomes 
invalid. Physically we might suspect that the bubble splits into two dis- 
joint parts. Second, the model can break down locally if a zero of d^,jf{uj,t) 
reaches the unit circle, which results in a cusp of the interface. It is well 
known that this is a common mechanism for breakdown in the unregularized 
model, (e = 0). 

Global breakdown has been observed with curvature regularization (see, 
e.g., |34j ) . and is also observed in our model. Figs.ll4landll5lshow examples, 
where each figure shows the evolution of a given initial condition for several 
values of e. For e = cusps do form, as is particularly obvious in Fig. [T5b . 
The regularization e > suppresses the cusps, but does not change the 
tendency to split into two parts. 

Whether local breakdown by cusp formation can occur in the regularized 
model, is a more difficult question. We recall that the neighborhood of 
a; = — 1 shows a special dynamics. The linearized evolution of the interface 
can lead to a very complicated shape near u; = — 1 since with increasing 
time singularities of /(a;,0) are gathered in this neighborhood. According 
to Subsect. I3DI these singularities for e > 1/2 dominate the local structure 
of the interface. We also note that for e = 00 even the linearized evolution 
is singular at a; = — 1, where it produces a spike. It thus is conceivable that 
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Figure 15: /(a;,0) = ^ + 0.75u; - 0.2u;=^ {N = 512, At = 0.0005). The 
solutions are shown at t = O.Oln for (a),(b),(c) and at t = 0.05 for (d). 

the nonlinear evolution yields a cusp or some other type of singularity at 
u = -1. 

We studied this problem with the initial condition /(w, 0) = l/u — O.lco'^. 
The results of the nonlinear evolution are shown in Fig. [161 Clearly the cusps 
forming for e = in the front part are suppressed for e > 0. We further 
observe that for e =1/100 or 1/10, the curvature near u = —1 decreases, 
whereas it increases for e = 1. This suggests that for e = 1 a cusp may be 
formed. 

5 Summary and conclusion 

Consistent with the eigenvalue analysis presented in [19], the results of the 
present paper strongly suggest that a uniformly translating circle is a lin- 
early stable solution of a Laplacian interface model regularized by a kinetic 
undercooling boundary condition. Furthermore, numerical results of the full 
nonlinear evolution indicate that the circle has a finite basin of attraction 
in a space of analytic functions. An important feature of the stabilizing 
mechanism is the advection of perturbations towards the back of the cir- 
cle. Except for a small region at the back that asymptotically contracts to 
a point, the final relaxation to the circle is exponential. With decreasing 
regularization parameter e > the anomalous behavior at the back is sup- 
pressed. However, perturbations increase as long as they are in the front 
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half of the circle, and this effect is strongly enhanced by lowering e. Since 
larger perturbations may lead to branching, this indicates that the basin of 
attraction of the circle shrinks exponentially with decreasing e. 

The interface model considered here is a reduced form of a PDE-model 
describing the streamer stage of electric breakdown in the simplest physi- 
cally relevant situation. It ignores the physics inside the streamer and the 
internal structure of the screening layer; the layer is approximated by the 
interface together with the boundary condition (l3|) which introduces the 
regularization. Numerical solutions of the PDE-model indicate that these 
approximations for sufficiently strong externally applied fields are justified in 
the dynamically active front part of the streamer. The back of the streamer 
is not represented adequately by the interface model. However, the evolution 
of the streamer and in particular stability or instability against branching is 
determined by the active head region, which corresponds to the front half of 
the circle in our analysis. Indeed, numerical solutions of the PDE-model in 
two dimensions show a behavior quite similar to the evolution of the front 
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half of weakly disturbed circles in the interface model. After reaching the 
streamer stage the streamer head is of nearly circular shape and moves with 
constant velocity. It slowly flattens at the tip and branches. Compared to 
the results of the interface model as illustrated in Fig. [T2l the main differ- 
ence is a slow increase of the head radius due to weak currents flowing into 
the head from the interior of the streamer. 

In summary, we believe that our results not only are of some interest in 
the context of interface models but also shed some light on the problem of 
streamer branching. 
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Appendix: Numerical calculation of the nonlinear 
evolution 

As explained in Sect. [2] the shape of the interface is given by 

z = X + iy = f (e*°, t) , — vr < a < vr. 

We restrict ourselves to interfaces symmetric with respect to the real axis, 
so that 

r (e*°,t) =/(e-^",t), 

with the corresponding equation holding for the potential $ (e*", t) . We use 
the Fourier representation 

oo 

f = Yl ^k^"", (93) 

oo 

^ = Cfc(i)e*'", (94) 

with a cutoff at high wave number k = N. Due to the symmetry, afc(t) and 
Cfc(t) are real, and the boundary condition at infinity ([5]) enforces c_i(t) = 
a-i{t). 

For a given shape of the interface the potential is determined by Eq. (jlOp : 
\daf {e'",t) I Re [$ (e*°, t)] = eRe (e'°, t)] (95) 

We represent \daf\ as 

oo 

\daf\ = dk{t)e'^'', (96) 

fc=— oo 

where the symmetry enforces = d^k £ I^- Foi^ a given / in Fourier 
representation (f93]) . 

oo 

5a/ = 5] iA:afc(t)e^'=°. (97) 
fc=-i 

The nonlinear term \dctf\ is computed via the standard pseudo-spectral ap- 
proach, i.e. \daf \ is obtained in the physical domain via inverse Fourier 
transform of Fourier coefficients in (j97p and taking the absolute value and 
then dk is determined by the Fourier transform of |5q,/|. Substituting 
Eqs. ([M|) . ([96]) into Eq. ([95]) . we find a system of linear equations for Cfc, 
A: > 0, which can be written as 

oo 

Y^idm-k + dm+k + ^m 5m,k)ck = {e5m,i-dm+i-dm-i)a-i, m>0. (98) 

fc=0 
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Here 6m,k denotes Kronnecker's symbol, and we used the identity c_i = a_i. 
We solve these equations with a cut off fc, m < N. Note that is needed 
up to A; = 2N. 

The evolution of the interface is determined by Eq. Q, which can be 
written as 



Re 



dtf 



Re [-ida^ (e*°)] 



R{a) 



(99) 



J is analytic for u € lA^ and is real for u; = by construction. Eq. (I99p 
therefore implies 



dtf{uJ,t) 



ujd^f{u;,t) 



1 

2^ 



■ UJ 



R{a') da', 



which for uj — > e^^ reduces to 



dtf 



ujdujf 



i a' , 

= R{a) P / cot — R[a + a) da' , 



(100) 



where P denotes the principle value. Symmetry enforces R{a) = R{—a), so 
that R{a) can be represented as 



R{a) = r,fc cos(A:a), G 



(101) 



fc=0 



where the again are determined by the Fourier-cosine transform numeri- 
cally. Substituting the expansions ([93]), (fTOT]) into Eq. (fTOOl) . we get 



dak 
dt 



fc+i 



^(fc - n)afc_„r„, k>-l. 



(102) 



n=0 



We again truncate this system of ODE's at A; = and solve it via 4-th 
order Runge-Kutta method (RK4) . Let the initial value problem (jl02p 
be specified as follows. 



^ = /(i,2/), y(io) = yo, 



(103) 



where y denotes the vector function (a_i, oq, ai, a^r). Then, the RK4 
method for this problem is given by the following equations 



Vn+i = yn + ^Hki + 2k2 + Sks + k^) 



tn+l = tn + h, 



(104) 
(105) 
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where yn+i is the RK4 approximation of y{tn+i), 



h = f{tn,yn), 

h = f{tn + -h, Vn + -hki) 

h = f{tn + ^h, Vn + ^hk2) 

h = f{tn + h,yn + hk3), 



(106) 
(107) 



(108) 
(109) 



and h is the time step. In the numerical implementation, h needs to be 
chosen small enough to ensure numerical stability and it is usually inverse 
proportional to the cut-ff A^. The cut-off needs to be chosen large enough 
so that the interface can smoothly represented, i.e. the Fourier coefficients 
are exponentially decayed for large k. For most of the initial conditions we 
used, there are only few Fourier coefficients are not zero. As time evolved, 
number of nonzero Fourier coefficients will increase. When the high fre- 
quency mode is no longer exponentially small, the algorithm needs to be 
terminated or more Fourier modes need to be used. Adaptive Fourier mode 
is beyond the scope of this paper. Here we only used fix cut-off N and make 
sure that the high frequency modes are exponentially small at later time. In 
the numerical simulations, we use both double and quadruple precision to 
compute solutions for large enough N. To prevent the spurious growth of 
the high-wavenumber coefficient generated by run-off error, we filter out the 
coefficient which is below the chosen threshold. If the threshold is chosen 
to be too large, aliasing may occurs. If the threshold is chosen to be too 
small, it cannot effectively reduce the run-off error. The reasonable choice 
from experience is about 1000 bigger than the run-off error. We choose the 
threshold to be 10~^^ for double precision and 10~^^ for the quadruple pre- 
cision. We can compare the results from both precision to ensure the results 
we obtained are not spurious. 

Notice that the numerical simulation need to stop at some finite time 
because of singularity. When singularity is developed, the numerical results 
become unreliable. There is a way to test the accuracy of numerical results 
without knowing the exact solution. Suppose the numerical method is of 
p-th. order, we expect that 



y{nh) 



0{hY. 



(110) 




(111) 
(112) 



Thus the order p can be estimated by using the formula 




~ VAu 
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We choose the initial condition f{uj,0) = — O.lu;^ and compute solutions 
for step size 0.001, 0.0005, and 0.00025 with N = 64. In Fig. [HI we can see 
that the order stays close to 4 up to the time equals to 2.5, 0.7, 1.25, and 0.9 
for (a) e=0 (b) e=0.01 (c) e=0.1 and (d) e = 1 respectively. Another way to 
test the accuracy is to check whether the area conservation holds or break 
down. The area enclosed by the interface should remain as a constant and 
it can be estimated by 

N 

A = -Y, kakitf. 

k=-l 

In Fig. [TSl the area changes versus time are shown for step size 0.00025. It 
is clearly that the area conservation and order of accuracy break down at the 
similar time for (a) e=0, (b) e=0.01, and (c) e=0.1. Notice that, in accuracy 
test, the order drops to the first order at t ~ 1 and the area conservation 
still holds up to i ~ 1.3. Similar behaviors are observed for other initial 
conditions and different Fourier modes A'^. For the computational results 
shown in the manuscipt, we show the sulotions up to the time that acuuracy 
of the solutions can be assured. 




0.5 1 1.5 0.5 1 1.5 



Figure 17: Numerical estimation of the order of accuracy for the initial 
condition /(u;,0) = ^ - O.lw^ with (a) e=0, (b) e=0.01, (c) e=0.1, and (d) 
e=l. (AT = 64) 
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0.5 1 1.5 



Figure 18: Estimation of area for the initial condition f{u!,0) = — — O.lo;^ 
with (a) e=0, (b) e=0.01, (c) e=0.1, and (d) e=l. {N = 64) 
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